High-spin binary black hole mergers 
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We study identical mass black hole binaries with spins perpendicular to the binary's orbital plane. 



These binaries have individual spins ranging from s/m 



-0.90 to 0.90, (si 



S2 in all cases) 



which is near the limit possible with standard Bowen-York puncture initial data. The extreme 
cases correspond to the largest initial spin simulations to date. Our results expand the parameter 
space covered by Rezzolla et al. and, when combining both data sets, we obtain estimations for the 
minimum and maximum values for the intrinsic angular momenta of the remnant of binary black 
hole mergers of J/M 2 = 0.341 ± 0.004 and 0.951 ± 0.004 respectively. Note, however, that these 
values are reached through extrapolation to the singular cases sij — j S2 1 
estimates until full-fledged numerical simulations provide confirmation. 



1 and thus remain as 



PACS numbers: 04.25.Dm, 04.70.Bw, 95.30.Sf, 97.60. Lf 



I. INTRODUCTION 

The existence of black holes, originally introduced as 
a family of solutions to the vacuum Einstein field equa- 
tions, was a matter of speculation for the best part of the 
20th century. In the past decade, however, astronomical 
observations placed them as the most promising models 
for objects detected in X-ray binaries (with sizes of a few 
to tens of solar masses) and for the supermassive entities 
at the center of some galaxies (with millions to billions 
of solar masses) [l|. Any black hole can be fully speci- 
fied by its mass, angular momentum and charge. Since 
electrically charged black holes are quickly neutralized 
by free charges found in their vicinity (i.e, from accre- 
tion disks, interstellar plasma, etc.), only their mass and 
angular momentum are of astrophysical relevance. In 
a dynamical environment, several factors can determine 
the black hole's angular momentum: the characteristics 
of its progenitor (in case of stellar collapse formation), 
the merger with other black holes and neutron stars of 
comparable mass, the merger with smaller objects, and 
accretion of matter from a surrounding disk. (For a re- 
view on how these situations can alter the black hole 
rate of rotation and the bounds for the maximum spin 
attainable see, for instance, and references therein.) 
Highly spinning black hole binaries are of particular in- 
terest given their astrophysical relevance. X-ray spectro- 
scopic studies of accretion disks around the supermassive 
galactic black holes may provide evidence of spin param- 
eters larger than 0.9 [|[ (however, see for instance [1] on 
how this may only be a tentative estimate) . 

Recent progress in numerical relativistic simulations 
of binary black holes (BBH) @,i,0] makes now possible 
long and stable evolutions that were impractical a few 
years ago. In this paper, we studied identical black holes 
with spins perpendicular to the orbital plane. In general, 
black hole spins would have arbitrary directions. How- 
ever, it has been recently suggested [8] that supermassive 
galactic BBH may favor configurations with spin align- 
ments like the ones studied here, due to the dynamical 



interaction of the holes with their galactic environments. 
Several groups have recently studied BBH with spins per- 
pendicular to the orbital plane [9j, [l0|, llll LL2|, |l3|. Here, 
we extend those studies, with binaries with initial indi- 
vidual spins ranging from s/m 2 = —0.90 to 0.90 (these 
extreme values are the largest modeled to date) with 
mi = m% and s\ = S2- The binaries in this sequence 
have not been considered before in such detail. Note also 
that, due to the symmetry of these systems, no gravita- 
tional recoil is present in the merger remnant. 

We perform a least-square data fit of our data following 
the quadratic formula used by Campanelli et al. 0] and 
Rezzolla et al. [l3[ and also of both Rezzolla et al. and our 
data sets combined. The last fit predicts minimum and 
maximum values for the spin parameter of the black hole 
remnant of 0.341 ± 0.004 and 0.951 ± 0.004 respectively. 

We find that current numerical techniques for the evo- 
lution of black holes with spins s/m 2 > 0.75 are limited 
in that they produce an artificial loss of angular momen- 
tum that increases with the magnitude of the spin. While 
this effect can be diminished by increasing the grid res- 
olution, even relatively large resolutions such as M/85 
present loss rates larger than 1% per 100M for the case 
s/m 2 — 0.90. This effect is clearly seen in the long term 
evolution of BBH that results in a highly spinning black 
hole and also on single black hole simulations. 

In Sec. [TTJ, we present the numerical details and tests of 
our simulations. Sections IIIII and [TVl present our results 
and conclusions. 



II. NUMERICAL TECHNIQUES AND TESTS 

All the binary systems studied here consist of identical 
black holes with spins aligned with the orbital angular 
momentum. In our coordinates, the orbit develops in the 
xy plane and the only non-vanishing component of the 
spins is in the z direction. We simulated binaries with 
spin parameters ranging from s/m 2 = —0.90 to 0.90. 

As starting points of our simulations, we use stan- 
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TABLE I: Initial data parameters. Here mj is the bare 
mass parameter of each puncture and M = 2m is the sum 
of the ADM masses m measured at each puncture. The 
holes have coordinate separation D, with puncture locations 
(0, ±D/2, 0), linear momenta (^fP, 0, 0), and spins (0,0, s). 
We also list the initial values of the ADM mass Mi and the 
angular momentum Ji. The 2PN angular velocity is set to 
SIM = 0.05550 in each case. 



dard puncture initial data [14| with the momentum 
and spin parameters in the extrinsic curvature given 
by second-order post-Newtonian (2PN) estimates [loT ]. 
It is sufficient to use 2PN estimates because standard 
puncture data are inconsistent with PN theory beyond 
(v/cf [H [13, El El. These parameters along with 
the initial ADM mass Mj and angular momentum Ji are 
shown in Table |TJ The coordinate distance D, the linear 
momenta P and spin parameters s are directly used in 
the Bowen-York extrinsic curvature, while the bare mass 
parameter is obtained from the condition that the ADM 
masses measured at each puncture should be m = M/2. 
As in [H Hi], H3| we assume that m is a good approxi- 
mation for the initial individual black hole masses. 

To complete the definition of the initial data, we also 
need to specify initial values for the lapse a and shift 
vector f3 l . At time £ = we use 

/ m b m h \~ 4 

01 = 1 + 1 1" ~A 

V 4n 4r 2 y 
P = 0, 

where ta is the coordinate distance from puncture A. 
Lapse and shift evolve according to 

(d t -(3 l di)a = -2aK, 
{d t -f3 k d k )p = \b\ 
{d t -(3 k d k )B l = (d t - (3 k d k )V - V B\ 

The gravitational fields are evolved using the "moving 
punctures" method [f| 0] with the implementation dis- 
cussed in [H, [24| , with the exception that as in [25j the 
dynamical variable <fi has been replaced by the variable 

W = e" 2 ^, 



\jT%Ob 


Structure 


1 


[5 x 48 


5 X 54] [hmin 


= M/56.9 


OB = 


238. 5M] 


2 


[5 x 56 


5 x 63] [h min 


= M/66.4 


OB = 


235.3M] 


6 


[5 x 64 


5 X 72] [h m in 


= M/75.9 


OB = 


246.4M] 


4 


[5 x 72 


5 X 81] [hmin 


= M/85.3 


OB = 


243.0M] 


5 


[5 x 80 


5 X 90] [hmin 


= M/94.8 


OB = 


240. 3M] 


6 


[5 x 80 


5 X 90] [hmin 


= M/85.3 


OB = 


267.0M] 


L 


[3 x 69 : 


6 x 149] [h min 


= M/80.0 


: OB = 


= 256. 0M] 



TABLE II: Characteristics of our numerical grids. The values 
between brackets show the number of inner (moving) levels 
times the number of grid points per level per dimension, plus 
the number of outer (fixed) levels times the number of grid 
points, hmin is the finest grid spacing and OB gives the coor- 
dinate distance to the grid's outer boundary. The last entry 
corresponds to a grid used with the code LEAN [26| . 



which obeys the evolution equation 

d t W = \w (aK - d. t f3 l ) + (3 i d l W. 

The new variable W has two advantages. First, our sim- 
ulations seem to converge better when we use W instead 
of (f> or x = e~ 4( ^ [6]. This may be related to the fact 
that W grows linearly near the black hole centers after 
some time of evolving the system. Such linear behav- 
ior leads to more accurate finite differencing derivatives 
near the punctures. In addition, the determinant of the 
3-metric det(7y ) = W~ 6 remains always positive, even 
if W becomes slightly negative due to numerical error. 
This property is not ensured if one evolves, for example, 
the variable \- I n the latter case one has to explicitly 
guard against this to prevent code crashes. 

The simulations presented here were performed using 
the BAM code, details of which are described in [23112^]. 
BAM is based on three-dimensional Cartesian coordi- 
nates and can achieve high spatial resolution near the 
punctures using a structure of nested refinement levels. 
The outermost of these levels are fixed, while the inner- 
most track the movement of the punctures. For the runs 
presented here, we used 10 levels of refinement with the 
outer boundaries located about 240M away from the grid 
center. Since the black holes are identical, we use quad- 
rant symmetry. Table |TT] lists the characteristics of the 
different grid layouts used in this paper. BAM character- 
istics make it ideal for quick and relatively inexpensive 
runs [24J. Most of the simulations presented here were 
done on dual processor workstations with characteristics 
speeds of 0.9, 1.7 and 2.5 days per orbit when using grids 
1, 2 and 3 respectively. 

Table [nil enumerates the binary simulations using the 
nomenclature SXX_Y, where XX indicates the value of 
s /m 2 and Y indicates which one of the numerical grids of 
TableHIlwas used. We measure the mass and angular mo- 
mentum of the remnant of the binary black hole merger 
using techniques based on volume integrals (see [24[ for 
details). These values of mass and angular momentum 
are labeled Mf and J/ respectively. 



3 





MtIM 

xvx j i xvx 


It IM 2 

Jf/lVl 


J i l al .f 


c n nn i 

b — u.yu_i 


u.y (o 


n Q A o 

0.348 


0.360 


5 — 0.90_2 


0.970 


0.359 


0.382 


c n nn q 


n 07n 
u.y * u 


U.ooo 


U.ooU 


b — U.9U_4 


0.970 


0.358 


0.380 


j — u. ( o_ i 


n 07n 
u.y / u 


U.4U0 


U.4ol 


5 1 — 0.75_2 


0.968 


0.408 


0.436 


S — 0.75_3 


0.968 


0.409 


0.437 


S — 0.5CL1 


0.963 


0.483 


0.520 


o — U.oU_z 


n n^Q 


U.486 


U.024 


S — 0.25_1 


0.958 


0.551 


0.600 


q n ok o 
o — U.zo_z 


n qkq 
u.yoo 


U.ooo 


n fine 


b + O.OtLl 


0.951 


0.605 


0.669 


c _i_ n nn o 
o + U.UU_z 


u.yoi 


U.014 


U.d / y 


c i n nn o 


n nc i 

u.yoi 


0.619 


n co a 
0.684 


c i n nn /i 

o + U.UU_4 


0.951 


0.619 


0.684 


c i n ok i 


n n/i a 
U.y 44 


n c*7i 
U.6 1 1 


U. <03 


c* i n o c n 

o + 0.25-2 


0.944 


0.676 


0.759 


S + 0.50_1 


0.934 


0.721 


0.827 


S + 0.50_2 


0.933 


0.720 


0.826 


c i n co i 
b + U.OZ_l 


0.92/ 


n to 1 
0. Ml 


(1 ori 

O.ool 


A + 0.62-2 


0.926 


0.737 


0.859 


A + 0.62_3 


0.926 


0.737 


0.859 


o + U. (0-1 


U.91 ( 


U. ( 29 


U.866 


A + 0.75-2 


0.916 


0.733 


0.874 


c~< i r\ Tr" o 

b + 0.75_3 


0.916 


0.736 


0.877 


C~< i r\ TV A 

b + 0.75-4 


0.916 


0.738 


0.880 


5 + 0.82-1 


0.911 


0.688 


0.828 


5 + 0.82_2 


0.910 


0.703 


0.850 


<? -|_ D £9 3 

> J U.OZj— o 


n QOQ 


n 7m 


fl 848 

U.O'iO 


S + 0.90.1 


0.905 


0.645 


0.788 


5 + 0.90_2 


0.905 


0.640 


0.781 


5 + 0.90.3 


0.903 


0.643 


0.788 


5 + 0.90_4 


0.902 


0.643 


0.790 


5 + 0.90.5 


0.902 


0.646 


0.794 


5 + 0.90.6 


0.906 


0.644 


0.785 



TABLE III: Binary black hole simulations performed for this 
paper. The runs are named SXXJY, where XX indicates 
the value of s/m 2 and Y is the numerical grid used for that 
run. Table ITT1 describes the different grids. 



We study the effect of different approximations and 
limitations inherent to our simulations: grid resolution, 
grid structure, measurement of the final mass and an- 
gular momentum, and the intrinsic characteristics of our 
initial data sets. We start by varying the spatial resolu- 
tion of our grids, while keeping the grid size and struc- 
ture (i.e., the layout of the refinement levels) unchanged. 
Grids 1 to 5 cover maximum spatial resolutions going 
from M/56.9 to M/94.8. The results, presented in Ta- 
ble 012 do not seem to vary significantly with the spatial 
resolution: the values of Jf/Mf obtained for grids 2 to 
5 agree with each other at a level of about 1.5% or bet- 
ter. We test the convergence of the calculations of mass 
and angular momentum, using four runs for the binary 
with spin s/m 2 = 0.90 (labeled S + 0.90.F, Y=l to 4). 
These runs shared the same grid layout but with vary- 
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FIG. 1: Mass and angular momentum plots for three different 
resolutions. The dashed curves have been scaled according to 
a factor corresponding to 4th order convergence in spatial 
resolution. 



ing spatial resolution, ranging from h min = M/56.9 to 
hmin — M/85.3 [36|] . We find that the lowest resolution 
run (5 + 0.90.1) fell outside the convergence regime. The 
results of the other three runs, presented in Fig. [1] seem 
to show that these runs are in the convergent regime. 

To evaluate the influence of the grid layout, we com- 
pare the results from two grids with identical resolution 
but different number of grid points per box (Grids 4 and 
6 from Table HH). The box size has been found crucial in 
previous work [23j : boxes smaller than critical sizes tend 
to change drastically the binary dynamics, altering orbits 
and merger times. Runs 5 + 0.90-4 and 5 + 0.90-6 re- 
turned values of mass and angular momentum with differ- 
ences of AM f = 4 10~ 3 (0.4%) and AJ f = 5 10- 4 (0.1%) 
respectively. 

Additionally, and as an independent check, we compare 
our results with those of Campanelli et al. [28J fshown in 
Fig. [2] as empty circles) and Pollney et al. [12[ (empty 
square). 

To evaluate the accuracy of the algorithms used to 
measure mass and angular momentum, we compare our 
volume-integral based results with calculations done us- 
ing surface integrals [24[ and find differences of up to 
0.5% in magnitude. We also test the satisfaction of the 
Christodoulou formula J c = 2 M irr (Mf-M^) 1 / 2 (M irr 
being the irreducible mass of the final black hole) . The 
relative difference in the values of J/ and J c for the run 
5 + 0.90_5 was less than 0.5%. 



III. RESULTS 

Figure [5] summarizes the results of Table IIIII The spin 
parameter of the binary remnant Jf/M 2 is shown as a 

function of the initial black hole spins s/m 2 . These val- 
ues were measured at 500M after the merger (which cor- 
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FIG. 2: Spin of the merger remnant as a function of the initial 
black hole spins measured 500Af after the black hole merger. 
The binaries are composed of identical black holes with ini- 
tial spins of magnitude s/m 2 aligned with the orbital angular. 
The values calculated in [2f| ([l^]) are shown as empty cir- 
cles (square). The solid square represents the remnant spin 
calculated right after the merger for the case s/m 2 = 0.90, 
using Grid L. The solid line is a quadratic interpolation of 
the values corresponding to Grid 2 for s/m 2 < 0.75 plus the 
Grid L value. 



responds approximately to the side-to-side light-crossing 
time of our grids). The most striking feature of the plot 
is the existence of an apparent maximum at s/m 2 ss 0.75. 
As we will see, this maximum is merely a numerical arti- 
fact caused by an artificial loss of angular momentum in 
the case of highly spinning black holes if we wait for 500AZ 
after the merger. Shortly after the merger, the spin of 
the remaining black hole is still larger. This is confirmed 
by a simulation of the s/m 2 — 0.90 case using the LEAN 
code (Grid L in Table [TTJ) that tracks the common ap- 
parent horizon and the emission of energy and angular 
momentum in the form of gravitational waves. We calcu- 
late the spin about 50M after the merger in four different 
ways: from the balance of gravitational wave energy and 
angular momentum loss, from the ringdown frequencies 
(using the tabulated results given in [29J]), from the iso- 
lated horizon geometry (using the techniques from 0]), 
and using Christodoulou's formula. The corresponding- 
results are Jf/Mj = 0.95, 0.925, 0.913, 0.918 respec- 
tively. The first result is the least accurate, given that 
the wave extraction radius (50M) was not considerably 
far from the center. The rest cluster around 0.92 (plotted 
in Fig. [2] as a solid square) that agrees better with an 
extrapolation from lower s/m 2 results. Figure [3] shows 
the evolution of the angular momentum obtained using 
isolated horizon techniques which indicates a loss rate of 
about 0.5% per 100M. This last result together with the 
measurements at 500M after the merger seem to indi- 



FIG. 3: Evolution of the angular momentum for the binary 
with s/m 2 — 0.90 obtained using the isolated horizon tech- 
niques of 0. The top plot corresponds to the angular mo- 
mentum of the remnant of the merger, while the bottom plot 
tracks the spin of the individual black holes before the merger 
(occurring at 265M). 



cate an artificial loss in angular momentum for high spin 
black holes. Also note that this effect is not seen at the 
lower spin end of Fig. O This is simply due to the fact 
that the merger remnant of binaries with individual spins 
s/m 2 < 0.75 is a low spinning black hole, and thus not 
affected by this effect. 

To study this effect in more detail, we perform single 
black hole evolutions. Figure 2] compares evolutions for 
the cases with s/m 2 = 0.53 and 0.90. In order to facili- 
tate the comparison, we plot the relative differences from 
the initial values of the mass, angular momentum and 
intrinsic angular momentum parameter, these quantities 
measured using volume integrals methods. The angular 
momentum loss, while negligible in the 0.53 case, is more 
pronounced in the larger spin case. Note however that 
the mass is largely unaffected in both cases. Figure [5] 
shows the results for four different grid resolutions, in- 
dicating firstly that the mass is better conserved than 
the angular momentum and secondly that this loss gets 
smaller with increasing resolution, albeit quite slowly. 
The angular momentum curves are still consistent with 
4th order convergence, however, Richardson extrapolat- 
ing these curves to arbitraily large resolution still shows 
loss of angular momentum. For BBH evolutions with 
many orbits before the merger, this loss could in princi- 
ple produce an "effective" value of s/m 2 at the time of 
the merger smaller than the initially intended. In order 
to evaluate this effect, we track the evolution of the in- 
dividual spins before the merger as shown in the bottom 
panel of Fig. [3] An angular momentum loss of about 3% 
is detected before the merger, which occurs at t = 2G5i\/. 




t/M 



FIG. 4: Single black hole evolutions with s/m 2 = 0.53 
(dashed) and 0.90 (solid), on Grid 4. The plots present the 
relative differences with respect to the initial values, denoted 
with the sub-index. The bottom plot presents the relative 
variation of intrinsic angular momentum j = Jf/Mr. 



FIG. 6: Puncture tracks of one of the holes for two simula- 
tions with s/m 2 — 0.90 that start out at different coordinate 
separations; D = 8M (dashed) and D = 6M (solid). The 
D = 6M orbits have been rotated to coincide with the ones 
from the D = 8M case in the final orbits. Both runs were 
performed using Grid 3. 
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FIG. 5: Evolution of single black holes with s/m 2 =0.90 for 
different spatial resolutions. 



At first sight, part of this loss could also be attributed 
to our choice of initial data. Traditionally, puncture ini- 
tial data sets are constructed by solving the Hamiltonian 
constraint under the assumption of conformal flatness for 
the spatial metric. The momentum constraint is analyti- 
cally satisfied by the use of the Bowen-York (BY) formula 
for the extrinsic curvature which, in principle, allows for 
the arbitrary specification of the linear and angular mo- 
mentum of the black holes. In practice, the amount of 



"junk radiation" associated with these initial data sets 
increases with the magnitude of the spin of the hole. The 
amount of this extra radiation has been studied pertur- 
batively by Gleiser et al. [30| and numerically by Hannam 
et al. [3l|] and Dain et al. [32]. The last shows that, for 
the case s/m 2 = 0.90, a single BY black hole would relax 
into a Kerr black hole after emitting about 0.1% of the 
initial gravitational mass. 

To further study this effect in the evolution prior to 
the merger, we perform a simulation of a binary with 
s/m 2 = 0.90 that starts out at a coordinate separation 
D = 8M and compare it with one starting at D = 6M. 
The former simulation lasted more than two orbits and 
a half longer than the latter, corresponding to an addi- 
tional simulation time of 395M. Figure[S]shows the orbits 
of one of the punctures for these simulations. Both runs 
were performed using Grid 3. Figure [7] shows the corre- 
sponding evolution of the mass and angular momentum, 
where the curves corresponding to the D = 8M run were 
shifted in time by 395 M . Naively, one would expect the 
D = 8M to dissipate more angular momentum before 
the merger, given that it is more than two orbits longer 
than its D = 6M counterpart. This would result in the 
longer run remnant having a smaller intrinsic spin than 
the shorter evolution. However, the resulting spin for the 
D = 8M run is about 3% larger for the D = 6M case, but 
within the error margin for these measurements, making 
it difficult to draw any conclusions. 

Next, we fit the values of the highest resolution runs 
from Table ITTT| ignoring the results for s/m 2 > 0.75, but 
adding the result at s/m 2 = 0.90 from Grid L (solid 
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FIG. 7: Evolution of mass and angular momentum for the 
runs of Fig. [6] The solid (dotted) lines corresponds to the 
separation D — 6AI (D — 8M). The dashed line corresponds 
to the D = 8M run shifted in time by 395M. 



remnant predicted by the fit of the combined data sets 
is J//M| = 0.341 ± 0.004 and 0.951 ± 0.004 respectively 
[3^ . Note that, while our data set is smaller than the 
one in [l3T |. it contains measurements that are closer to 
the extrapolated values for critical BBH. 

Recently, Boyle et al. [33[ introduced a formalism 
that predicts any final quantity resulting from the merger 
using a Taylor expansion on the initial binary mass ratio 
q = m\/m 2 and the components of the initial spins. For 
the highly symmetric binaries considered in this paper, 
their expansion for Jf /M 2 to second order reduces to the 

polynomial of Eq. ([T]) with the equivalences po 
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Pi 

P2 

Max J/M 2 
Min J/M 2 


0.6879 
0.1476 
-0.0093 
0.946 
0.355 


0.6883(4) 
0.1530(4) 
-0.0088(5) 
0.959(2) 
0.347(2) 


0.6855(16) 
0.1499(8) 
-0.0110(8) 
0.941(6) 
0.342(6) 


0.6888(4) 
0.1525(5) 
-0.0106(5) 
0.951(4) 
0.341(4) 



TABLE IV: Least-squares fit of Eq. (|TJ from Campanelli et 
al. (C), Rezzolla et al. [l3|] (R), this paper (M), and Rez- 
zolla et al. and our data sets combined (R+M). The last 
two rows show the extrapolation to initial critical black holes 
aligned (counter-aligned) with the orbital angular momentum 
that corresponds to the maximum (minimum) possible intrin- 
sic spin for the remnant black hole. 

square). We follow the fitting formula for the final black 
hole intrinsic spin used in Campanelli et al. Q and Rez- 
zolla et al. [Hi (Eq. (8)) 

Jf/M]=p + Pi {s 1 /m\ + s 2 /ml) 

+ P2 {s 1 /m 2 + s 2 /m 2 ) 2 , (1) 

where po. Pi and p 2 are the fitting parameters. We 
present in Table IIVI a comparison of our best fit param- 
eters with those from [9j and [l3j and with the ones ob- 
tained by fitting [l3| and our data sets together. For 
these fits, we assume a nominal error of 0.02 for all our 
values. We see that the fitting parameters are in an agree- 
ment consistent with small-number statistics in all cases. 
Table IIVI also shows the extrapolation to maximally ro- 
tating black holes, aligned and counter- aligned with the 
orbital angular momentum. The maximum and mini- 
mum values for the intrinsic angular momentum of the 



IV. CONCLUSIONS 

We studied the effect of the initial spins of the black 
holes in a binary system on the mass and angular momen- 
tum of the black hole that results from the merger. We 
concentrated on equal-mass binaries with spins aligned 
with the orbital angular momentum (s\ — s 2 ), covering a 
range of initial spin parameters going from s/m 2 — —0.90 
to 0.90. The runs at the extrema of the range are the 
highest spin simulations to date. The main results of the 
paper are presented in Fig. [21 where the spin parame- 
ter of the remnant (Jf/M 2 ) is given as a function of the 
initial spin parameters. 

We combined our results with those of Rezzolla et 
al. [iljj in a quadratic least-square fit and obtain, by ex- 
trapolation to initial critical black holes, predicted max- 
imum and minimum values of J/M 2 for the black hole 
remnant of 0.951 ± 0.004 and 0.341 ± 0.004 respectively. 
These error bounds are simply derived from the uncer- 
tainty of the fitting parameters provided by their stan- 
dard error. The small size of the samples studied here 
plus the fact that the limits to critical black hole results 
are obtained through extrapolation to points in parame- 
ter space with singular properties may lead to revisions 
in these estimates once binaries with spins larger than 
0.90 can be accurately simulated. For this, new recipes 
for initial data sets that allow for initial black holes with 
spin parameters larger than 0.928 (the limit of BY data) 
are needed. The methods introduced by Dain [3J, [35| 
and studied in head-on BBH simulations by Hannam et 
al. [3l| may hold the key to this problem. 

We also find a problem for the simulations starting 
with spins s/m 2 > 0.75. Current evolutions based on 
the moving punctures methods present losses of angular 
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momentum at non-negligible rates for highly spinning 
black holes, even when relatively high grid resolutions 
are employed. Measurements of the merger remnant in- 
trinsic spin 500M after the merger of two s/m 2 = 0.90 
black holes show the spurious loss of more than 10% of 
the angular momentum. This effect increases with the 
magnitude of the black hole spin. This loss does not 
affect strongly the calculation of gravitational wave tem- 
plates, since they cover only up to a short period after 
the merger. However, it becomes more important, for 
instance, for simulations of black hole-neutron star or 
neutron star-neutron star binaries, where the potential 
formation of an accretion disk around a black hole has to 
be followed for much longer periods. 
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